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Weakening dark-matter cusps by clumpy baryonic infall 
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ABSTRACT 

We consider the infall of a massive clump into a dark-matter halo as a simple and extreme 
model for the effect of baryonic physics (neglected in gravity-only simulations of large-scale 
structure formation) on the dark-matter We find that such an infalling clump is extremely 
efficient in altering the structure of the halo and reducing its central density: a clump of 1 % the 
mass of the halo can remove about twice its own mass from the inner halo and transform a cusp 
into a core or weaker cusp. If the clump is subsequently removed, mimicking a galactic wind, 
the central halo density is further reduced and the mass removed from the inner halo doubled. 
Lighter clumps are even more efficient: the ratio of removed mass to clump mass increases 
slightly towards smaller clump masses. This process is the more efficient the more radially 
anisotropic the initial dark-matter velocities. While such a clumpy infall may be somewhat 
unrealistic, it demonstrates that the baryons need to transfer only a small fraction of their 
initial energy to the dark matter via dynamical friction to explain the discrepancy between 
predicted dark-matter density profiles and those inferred from observations of dark-matter 
dominated galaxies. 

Key words: stellar dynamics - methods: A^-body simulations - galaxies: kinematics and 
dynamics - galaxies: structure - galaxies: haloes 



1 INTRODUCTION 

The success of the ACDM paradigm of cosmological struc- 
ture formation in reproducing the observed structure of the uni- 

■ verse on scales > 1 Mpc i s now well-established (see e.g. 

' ISpringel. Frenk & Whitell2006l and references therein). However, 
on the scale of individual galaxies (< lOOkpc) the density pro- 
files of dark matter haloes pose a potentially significant problem. 
The form of these profiles in the absence of baryonic physics 
has been extensively studied by means of numerical simu l ations . 
iDubinski & Carlberd ( Il99lh and iNavarro. Frenk & Whitel d 19971) 
showed that haloes on a range of mass scales have similar profiles, 
with the density in the innermost regions exhibiting an p oc r"' 
cusp. More recent work has found that the haloes can be better 
represe nted by profiles with either slightly shallower inn er cusp 
slopes jPehnen & McLaughlirjboOSl i lNavarro et al.ll2010h or con- 
tinuously varying slope, for example the lEinasto & Einastol ( 1 19721) 
profile (Navarro et al. 

However, there is mounting observational evidence that galax- 
ies occupy dark-matter haloes with almost uniform density cores 
with the stro ngest evidence coming from low surface brightness 
galaxies (see lde Bloldl200 9'. for a recent review). Although there 
have been fe wer such stud i es in h igh surface brightness spirals, re- 
cent work bv lSDanoetai] ( l2008l) has shown that in these galaxies 
also cored haloes are favoured. Even in low-luminosity, pressure 
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supported systems such as the dwarf spheroidal satellite galaxies 
surrounding the Milky Wa y, there is circumsta ntial evidence that 
their haloes are not cusped dOilmore et al.ll2007l) . 

It is important to remember however, that the ignorance of 
baryonic physics in the aforementioned simulations constitutes a 
significant limitation with regard to any discussion of the inner halo 
profiles. The inclusion of baryon physics is widely recognised as 
a crucial, albeit extremely technically challenging prerequisite for 
further progress in understanding galaxy formation. 

There are multiple ways in which the baryons can affect the 
density profile of a dark matter halo. First, a cloud of gas which 
initially extends throughout a cusped dark matter halo can dissi- 
pate energy by radiation and contract to the centre of the halo. As 
the dark matter responds to the deeper gravitational potential, this 
in turn leads to a ste epening of the dark m atter cusp, known as 
adiabatic contraction telumenthal et al.llT986) . This is not neces- 
sarily the end of the story, however, as baryons also have the ability 
to generate mass outflows driven by stellar winds and supemovae 
produced as a result of star formation. Depending on the efficiency 
of star formation, such processes can expel a significant fraction 
of the baryonic mass from the central regions of a galaxy, result- 
ing in a large-scale rearrangement of the d ark matter. Using semi- 
analytic arguments. 'Gnedi n & Zhao| ( l2002h found that when stellar 
mass loss was preceded by adiabatic contraction, the resulting halo 
distri bution was almost unc hanged from its original cusped pro- 
file. Read & Gilmord feOOSi) . however, subsequently showed that 
repeated episodes of adiabatic contraction followed by rapid mass 
expulsion could give rise to a reduction in the central density and 
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hence produce cored haloes. This occurs because although the ini- 
tial infall/outflow produces the same mass density profile, the ve- 
locity structure of the dark matter halo after the outflow is biased to- 
wards radial orbits. As a result, subsequent events are able to lower 
the inner density more easily. 

Another way in which the baryonic component of a galaxy 
can transfer energy to t he dar k matter halo is by means of a stellar 
bar. Weinberg & Katzl( l2002h proposed that the rotation of bars is 
decelerated by the exchange of energy and angular momentum with 
dark-matter particles o n orbits in reso n ance with th e bar's rotation. 
While further work bv ' Athanassouial ilOOA [2003b confirmed this 
conclusion, the size required for a stellar bar to significantly change 
the mass distribution in the inner regions of a dark-matter halo was 
found to be much larger than those observed in barred galaxies. 
Hence, this mechanism is thought to be of limited importance for 
the majority of dark-matter haloes, though the deceleration process 
is certainly affecting the evolution of galactic bars. 

However, galaxy formation also involves violent processes, 
where baryonic inflow is neither smooth nor adiabatic. Gas accre- 
tion is likely to occur during galaxy mergers, when it takes the form 
of clumpy infall rather than the slow contraction of a smooth cloud. 
Although a baryon clump falling into the centre of a dark halo will 
add to the gravitational potential there (and thus increase the bind- 
ing energy of the dark halo), it can, during this process, lose its or- 
bital e nergy via dynamical friction. lEl-Zant. Shlosman & Hoffmanl 
showed that the energy thus gained (i.e. binding energy lost) 
by the dark matter ca n produce an observ able impact on the dark 
matter density profile. lEl-Zant et al.l ( l2004l) and lNipoti et al.l ( l2004l) 
extended this work by performing A'-body simulations of galaxy 
clusters, where the infalling galaxies play the role of clumps. They 
found that the initial dark-matter cusp can be softened through the 
transfer of energy from the baryonic clumps to the dark-matter 
though the overall density profile, including the baryonic compo- 
nent, remained cusped. 

Attempts have also been made to add full baryonic physics 
to the studies of sinking clumps mainly based on the results of 
cosmological simulations, which try to model the e ffect of cool 



ing, metal enrichment and supernova feedback (e.g. [Gnedin et al 
2004; Romano-Diaz et al. 2008; Pedros a. Tissera & Scannapieco 
_2009; Johansson, Naab & Ostriker 2009). These simulations con- 
firmed the transfer of energy and angular momentum from baryons 
to the dark-matter, while the results on the dark-matter density re- 
duction were conflicting. This is presumably because of varying 
amounts, depending on the details of the respective model, of con- 
traction owed to the additional gravitational pull from the accreted 
baryons. 

Despite these promising attempts, a truly realistic modelling 
of baryonic physics is still beyond contemporary simulation tech- 
niques, not least because important baryonic processes, such as re- 
ionisation as well as primordial and ordinary star-formation, are 
themselves not sufficiently understood. Therefore, it is important 
to understand more quantitatively the purely stellar dynamical as- 
pect of this mechanism, which alone affects the dark-matter distri- 
bu tion. Some progres s towa rds t his goal has been m ade recently 
bv lJardel & Sellwood ( l2009h and lGoerdtet"ai1(l2010l) . However, a 
complete understanding of the pure stellar dynamical problem is 
still missing, but seems essential before attempting to interpret the 
results of simulations which include baiyons. In the present paper, 
we build on this previous work to explore the impact of clumpy 
baryonic infall more broadly. We consider more realistic initial con- 
ditions. In particular, we focus on clumps initially on parabolic or- 
bits, which may become bound to a halo during a merger event. 



and haloes with anisotropic velocity distributions, the expected sit- 
uation to obtain within the hierarchical structure formation scenario 
of CDM. 

The outline of the paper is as follows. In Section|2]we consider 
analytical estimates for the damage done to the halo by the energy 
transfer from the satellite orbit. Section |3]gives our modelling ap- 
proach for the A'-body simulations, while Sections|4]and|5]describe 
the resulting orbital decay and typical changes induced in the sim- 
ulated haloes. Sections |6] and |7] discuss the effect of halo velocity 
anisotropy on both orbital decay and damage to the halo, while Sec- 
tion [8] considers the effect of satellite mass and size. In Section |9] 
we demonstrate the effect of removing the accreted clump, corre- 
sponding to a galactic wind subsequent. Finally, Sections [TOland[m 
summarise and discuss our findings, respectively. 



2 THEORETICAL ARGUMENTS 



IChandrasekharl ' s (1943) dynamical friction formula for systems 
w ith a Maxwellian velocity distribution of dispersion a (eq. 8.7 
of lBinnev & Tremainell2008h 
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(In A is the Coulomb logarithm, and are the mass and velocity 
of the clump or satellite, and p is the mass density of dark-matter 
particles) shows that the deceleration is proportional to m^, such 
that the time for the orbit to decay /infaii nj^'. In particular, for 
this orbital decay to occur within (less than) a Hubble time, a mass 
of m, ~ 10*"** Mq, depending on the size of the dark-matter halo, is 
required. Chandrasekhar's formula also suggests that the drag force 
is strongest for small (because this increases the interaction time 
between perturber and dark-matter particles) and/or for large p. 

However, the formula cannot be used to assess the effect the 
infalling clump has on the dark matter. A simple estimate for 
the mass removed from the inner parts of the dark-matter halo 
can be obtained fr om the following argument originally due to 
dCoerdt et al.ll201ol preprint version). Assuming a circular orbit 
for the perturber, the specific energy lost when sinking from radius 
r -I- (5r to r is 
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with p(r) the mean density interior to radius r. Assuming this en- 
ergy is injected into the spherical shell at radius r, each dark-matter 
particle at that radius gains specific energy {Ss^l 8r)(m^ l Anpr^). A 
density core forms and the sinking of the clump stalls jRead et al.l 
l2006l) as soon as this energy equals the specific kinetic energy of 
each particle, which may be estimated as = GM(r)/2r. With 

(O, this yields (withy = -dlnp/dlnr) 



M(r) ~ 
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(3) 



This argument suggests (i) that the mass ejected by the perturber 
is comparable to its own mass, and (ii) that the density core which 
forms in response to the heating induced by the sinking baryonic 
clump has radius comparable to that at which the originally en- 
closed mass equals lUs. 

Strictly speaking, this argument only applies to circular orbits, 
which are not very realistic, and the assessment of the heating re- 
quired to turn a cusp into a core is rather crude. We now present 
a more quantitative estimate based on the exact energy difference 
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Figure 1. Predictions, using equation {S), for tlie ratio of the maximum excavated mass to tlie mass of the infalling clump (top), the radius at which this occurs 
(middle), and the radius inside of wliich half the mass has been removed (bottom; in units of halo scale radius) as function of the mass of the accreted clump. 
The clump was assumed to be initially either on a bound orbit with energy equal to that of the circular orbit at halo half-mass radius (left) or on a parabolic orbit 
(right). The initial halo is modelled to have density distribution )12t also used in our simulations, while for the final halo we used the models of equation (TT) 
(solid) or not (dashed: rj = 0.9 or dash-dotted: rj = 1.5). The different lines correspond to inner density slopes for the final halo between yo = (red) and 
yo = 0.4 (green), while along each line the scale radii ro of the final model are varied (and nij obtained from equation|8]. 



between initial and final halo and on the assumption that the orbital 
energy lost by the clump is absorbed by (the inner parts of) the halo. 

Assuming spherical symmetry, let M[(r) and Mf(r) denote 
the cumulative mass profiles of the initial and final halo. At large 
radii the halo is hardly altered, i.e. Mf(r) « M\(r), while at small 
radii the halo has been heated resulting in an expansion and hence 
Mf(r) < Mi(r). A quantitative relation between the change in M(r) 
and the mass of the clump can be obtained by considering the 
total energy budget. By virtue of the virial theorem, the total en- 
ergy of the initial equilibrium halo is half its potential energy Vj, 
to which the kinetic energy of the clump and the interaction energy 
between clump and halo must be added to obtain the total energy 
of the initial state (neglecting the clump self-energy) 



£i = i^(.? + cl.,(^,))4^.- 



(4) 



with r\ and the initial radius and speed of the clump and <I)i(r) 
the potential due to the initial halo. Expressing this in terms of the 
initial specific orbital energy Es = + of the clump and the 
cumulative mass profile gives 
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where we have used the relation V - 



^ M^(r)r-^ dr. For the 
final state (halo in equilibrium with the clump at rest in the centre) 
we have 
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- this is obtained as ii = V/2 using M(r) = + Mf(r) (assum- 
ing a point-mass clump) and ignoring the clump self-energy. If the 



clump is extended with cumulative mass profile ms(r), then <I>f(0) 
in equation ^ has to be replaced by 
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This equation relates the clump mass and its initial conditions with 
the change in the halo mass profile. Because this latter is a one- 
dimensional function, while relation ([8} provides only a single con- 
straint, it can be satisfied by many possible functional forms for the 
final mass profile Mf(r). However, using some simple yet reason- 
able models for the final mass profile Mf(r) we can obtain some 
quantitative estimates for the amount of mass excavated 



AM(r) = Mi(r) - Mf(r), 



(9) 



in particular its maximum and the radius at which it occurs, and 
their dependence on clump mass and initial specific energy e^. 

In order to compare directly to our simulations, we chose the 
same initial halo density profile, given in equation ( I12t below, as 
used in the simulations. For the final halo we assume two difi^erent 
families of models; the first have mass profile 



Mf(r) 



Mi(r) wilh = r" + r' 



(10) 



where yi is the central density slope of the initial halo. These mod- 
els have a p oc r^™ density cusp at small radii. The second family 
of models also have the central density slope yo and scale radius ro 
as free parameter and are given by 



Mf(r) = Mi(r)(l - exp(-[r/ro]'''-''»)). 



(11) 
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The top panels of Fig. [T] show the resuUing relations between 
clump mass rris, obtained for the above models from equation ((S), 
and the ratio of the maximum excavated mass over the clump mass, 
AM|^ax/™s for a clump on an orbit with initial specific energy of 
Es= - 0.12 (left, corresponding to a circular orbit at the halo half- 
mass radius), or an initially parabolic orbit (e^ = 0, right) decaying 
in a halo with initial density il2\ . as used in our A'-body simulations 
below. While these two types of models result in a range of values 
for the excavated mass, both models suggest that a clump of 1% of 
the halo mass can excavate about 1.7 times its own mass for e, = - 
0. 12 and twice its own mass for = 0. The models of equation Jilt 
also indicate that a less massive clump is relatively more efficient 
in that it removes a larger multiple of its own mass. 

The radius rn,ax at which AM(r) becomes maximal (and hence 
Pi = Pf ) is shown in the middle panels of Fig. [T] again as function 
of clump mass mj. For a clump of only 1% of the halo mass, this 
radius can easily reach the dark-matter scale radius. We also plot in 
the bottom panels the radius r5o% inside of which half the mass has 
been removed, i.e. Mf(r5o%) = ^Mi(r5o%), which is readily mea- 
sured from our A'-body simulations below and provides a more di- 
rect measure of the size of any possible density core. The rankings 
between the models of equation (TTJ {solid curves in Fig. [TJ w.r.t. 
''max and rjo* are opposite to each other: the model with the small- 
est r^o% has the largest r^^^ and 70 and vice versa. While both rso% 
and r,„ax clearly increase with ra^, as one would expect, this depen- 
dence is rather weak, approximately r oc mf^ for the models of 
equation dl It near = O.OlMhaio- 

Note that in obtaining these estimates we have assumed that 
the final halo is spherical and in perfect equilibrium. In reality, the 
momentum of the infalling clump is absorbed by the inner halo, 
which as a consequence moves slightly w.r.t. its outer parts. This 
motion is only weakly damped and it takes some time before its 
energy is transformed into internal heat. This implies that the es- 
timates for AMniax from equation ([8} (as used in Fig. [TJ are pos- 
sibly somewhat too high, depending on the amount of energy ab- 
sorbed into such oscillations. The inner asphericity resulting from 
the absorbed momentum also implies that the simple model spheri- 
cal over-estimates AMn,ax (because it under-estimates Ef given Mf). 



3 MODELLING APPROACH 

The above energy argument is only suggestive and cannot predict 
the final mass profiles and its dependence on the details of the 
satellite orbit and the initial halo equilibrium. To this end numeri- 
cal (A' -body) simulations are required. Recentlv. lJardel & Sellwooj 
1 I2OO9 ') performed such simulations for a satellite starting on a cir- 
cular orbit at the half-mass radius of a spheri cal dark-matter halo 
model with a Navarro, Frenk & White! d 19951) density profile and 
isotropic velocities. They found that the orbital decay of a satellite 
with 1% of the halo mass heats the dark matter and results in a den- 
sity core of radius ~ O.lro, where is the initial scale radius of the 
halo, the radius at which y (r) = -dl np/d In r = 2. 

A study by iGoerdt et al.l ( 120101) also considered satellites on 
circular orbits, but for a variety of halo density profiles with differ- 
ent inner density slope 7. These aut hors were especially interested 
in the stalling of the orbital decay d Ooerdt et al.ll2006l : iRead et al] 
[2OO6), which depends sensitively on 7. 

Assuming a circular orbit for the satellite may simplify some 
analytical arguments, but is certainly not very realistic. Similarly, 
the assumption of velocity isotropy for the dark matter is not jus- 
tified and made only for convenience (so that initial conditions are 



easily prepared), though simulations bv lArena & Bertir] (|2007|) in- 
dicate that velocity anisotropy plays no significant role. 

The aim of the present study is to extend the aforementioned 
simulations to more realistic initial conditions. In particular, we are 
interested in non-circular satellite orbits and cosmologically mo- 
tivated velocity anisotropy for the dark matter Our energy -based 
argument of ^ shows that the total amount of heating induced by 
the decaying satellite orbit only depends on the satellite orbit's en- 
ergy, but not on its eccentricity. However, where the satellite dumps 
its energy, and consequently which halo particles gain energy and 
angular momentum, will depend on eccentricity. In order to differ- 
entiate these effects, we will present two sets of A'-body simulations 
with satellite orbits of t he same energy. The first set essentially ex- 
tends the simulations of ljardel & Sellwood bv considering satellite 
orbits of varying eccentricity, but same initial radius and energy 
as used by those authors. The second set of simulations considers 
parabolic orbits, whose initial orbital energy just vanishes. If it were 
not for dynamical friction, these orbits would just pass by the dark- 
matter halo and never return. However, due to dynamical friction 
they become bound and decay to the halo centre if initially aimed 
sufficiently close. We consider isotropic velocities for the dark mat- 
ter as well as velocity anisotropy of various degree and radial de- 
pendence. Furthermore, we investigate the effects from changing 
the satellite mass and/or its adopted size. 



3.1 The halo model 

For the densit y profile of the dark-ma t ter ha lo, we adopt a trun- 
cated spherical lOehnen & McLaughlinI ( |2005|) model, which gives 
an excellent fit to simulated CDM halos and has density 

p(r) oc r-"" + sech(r/r,). (12) 

For Tt — > 00, this profile asymptotes to p oc r"''' and r"'''' at 
small and large radii, respectively, with a very smooth transition. 
We identify the scale radius with the radius at which -y(r) = 2 
(in the limit — > 00), which for these models is given by r2 = 
(11/13)''''' J ~ 0.687.S and set the truncation radius to rt = 10r2, 
which we identify with the virial radius. We consider various ve- 
locity anisotropy profiles; in particular models for which 

P-=^-^ (13) 

is constant and models for which 

4/9 

7777747?- (14) 

These latter models are isotropic in the centre and become increas- 
ingly radially anisotropic (for > 0) at large radii, again with 
a very smoot h transition, and are excellent descriptions of A'-body 
CDM haloes jPehnen & McLaughliij[200i) . 

To generate initial A'-body conditions for the halo, we sample 
positions from ( I12t and velocities from self-consistent distribution 
functions of the form L^^^f(E) for constant models with /(e) ob- 
tained from an Abel inversion (Cuddeford 1991). For models with 
j8(r) as in equation l ll4t . we generate initial conditions using the 
made-to-measure A'-body method of iDehnenI 12009). 

For models with constant yS, the resolution in the inner parts is 
enhanced by increasing the sampling probability by a factor g(e)"' 
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which is compensated by setting particle masses /j, proportional to 
giSi). We used 



g(e) 



1 

-t circ ' 

circ ^ ■' 



(15) 



with q=lO the ratio between maximum and minimum particle 
mass and rri,c(e) the radius of the circular orbit with specific en- 
ergy E. The gravitational forces were computed using a softening 
kernel with density profile given in equation ( 116) below and re- 
placed by the softening length e = 0.005. Testing this method for 
our particular purposes we found that it allows a reduction of N to 
half at the same central resolution without any adverse effects. 

We use a unit system where G = M = r2 = I, which implies a 
time unit of 15.5 Myv{r2/500pcf'^{M/lO^Me)-^'^. 



3.2 Orbital and other parameters of the infalling clump 

Henceforth, we shall use the term 'satellite' for the infalling bary- 
onic clump, which we model as a single massive extended (soft- 
ened) particle. For the satellite mass iris we considered 0.3%, 1%, 
and 3% of the total halo mass (only 77% of which is inside rt as- 
sociated with the virial radius), while the satellite size rs was taken 
to be 0.01, 0.03, or 0.1 times the halo scale radius. This means that 
the satellite is efi^ectively modelled to have spherical density profile 



15 
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8/r (r2 + 7-2)7/2 • 



(16) 



We considered a large range of initial satellite orbits, but report 
here only on two sets of simulations. The simulations of the first 
set, summarised in the top half of table [T] all start from r = 4, the 
radius containing 40% of the total or 54% of the mass within rt, and 
have specific orbital energy = -0.12 equal to that of the circu- 
lar orbit at that radius. The only remaining free parameter of these 
orbits is the pericentric radius of the initial orbit — owing to dynam- 
ical friction, the actual trajectory of the satellite may have a slightly 
smaller firs t pericentric radius. These simulations thus extend those 
reported bv lJardel & SellwoodI (l2009t) to non-zero eccentricity and 
also anisotropic halo velocity distributions. The second set of simu- 
lations, summarised in the second half of table[T] employs parabolic 
orbits, i.e. with initial orbital energy Es = 0, started at r = r,, cor- 
responding to the halo virial radius. Again, the only free orbital 
parameter is the pericentric radius of the initial orbit. 

Within either set of simulations the specific energy of the ini- 
tial orbit is the same. This choice was motivated by the analytic 
argument of section |2] which suggested that orbital energy is the 
main parameter afi'ecting the amount of 'damage' done to the halo. 
Thus keeping this energy fixed and varying orbital eccentricity or, 
equivalently, the pericentric radius, we can study the influence of 
this secondary orbital parameter 



Table 1. Initial conditions and results for our simulations. Initial conditions 
are specified by the satellite size Ts and mass ffis, which default to Ts = 0.03 
and nis = 0.01, respectively; the peri-centric radius rp^jj of the initial satellite 
orbit; and the halo initial velocity anisotropy /Jj , which is either constant or 
/3(r) given by equation U4t with /3„ = 1 . As results we list the time fi„faii 
for the satellite to fall to the centre of the halo (defined in Section 14. U ; 
the radius r^ax of maximum halo-mass reduction AM„,ax (compared to the 
control simulation); the radius rso% where the cumulative mass is reduced 
to 50% compared to the control simulation; the radius rM=m^ within which 
the final halo mass equals His; and the maximum (over all radii) of p/cr^ for 
the final halo. 



simulations started at = 4 and run until t = 1000 
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0.3 


145.9 


0.641 


0.0088 


0.127 


0.262 


2.46 


default 


4.01 


0.3 


178.8 


0.665 


0.0082 


0.128 


0.257 


2.37 


rs =0.01 


4.01 





181.8 


0.802 


0.0084 


0.048 


0.247 


4.67 


rs =0.1 


4.01 





231.1 


0.464 


0.0076 


0.215 


0.272 


2.60 


ms=0.003 


4.01 





491.4 


1.932 


0.0028 


0.042 


0.116 


13.38 


ms=0.03 


4.01 





90.3 


1.278 


0.0210 


0.257 


0.541 


0.88 




simulations started at 


ri = 10 and run until f 


= 2000 




rs, »is 


rperi 


yS, 


finfall 


' max 


AMn,ax 


r50% 


^ M=ms 


max (5) 


default 


0.0 


-0.43 


225 


0.802 


0.0243 


0.530 


0.384 


0.86 


default 


0.4 


-0.43 


460 


0.891 


0.0208 


0.393 


0.332 


1.26 


default 


0.8 


-0.43 


835 


1.271 


0.0180 


0.207 


0.280 


2.96 


default 


1.3 


-0.43 


1365 


1.922 


0.0163 


0.094 


0.256 


4.32 


default 


0.0 





215 


0.701 


0.0229 


0.515 


0.378 


0.69 


default 


0.4 





408 


0.734 


0.0197 


0.447 


0.352 


0.82 


default 


0.8 





774 


0.906 


0.0167 


0.331 


0.315 


1.20 


default 


1.3 





1319 


1.288 


0.0138 


0.265 


0.294 


1.62 


default 


0.0 


0.3 


202 


0.710 


0.0220 


0.511 


0.380 


0.57 


default 


0.4 


0.3 


388 


0.721 


0.0193 


0.445 


0.355 


0.69 


default 


0.8 


0.3 


700 


0.750 


0.0178 


0.389 


0.339 


0.88 


default 


1.3 


0.3 


1165 


0.853 


0.0160 


0.358 


0.329 


0.84 


default 


0.0 




183 


0.858 


0.0187 


0.445 


0.372 


0.52 


default 


0.4 


m 


320 


0.812 


0.0182 


0.435 


0.368 


0.53 


default 


0.8 


m 


593 


0.807 


0.0157 


0.379 


0.350 


0.66 


default 


1.3 


m 


947 


0.844 


0.0159 


0.372 


0.348 


0.65 


rs =0.01 


0.4 




249 


0.834 


0.0174 


0.412 


0.360 


0.58 


rs =0.1 


0.4 


m 


566 


0.822 


0.0177 


0.415 


0.360 


0.65 


ms=0.003 


0.4 




1902 


0.406 


0.0064 


0.207 


0.180 


2.02 


ms=0.03 


0.4 


m 


95 


1.149 


0.0427 


0.871 


0.762 


0.16 



3.3 Technicalities 

The A^-body simulations are performed using the public A^-body 
code gyrfalcON which uses the 0(N) force solver falcON jPehnenl 
l2002j) with minimum opening parameter 6„^;„ = 0.5 and employs 
an adaptive time-stepping scheme. 

The simulations were run for 1000 or 2000 time units for the 
first and second set of simulations, respectively. The energy con- 
servation was typically 1 part in 10* (more accurate control sim- 
ulations obtained the same results). Halo models with constant 



fi had 1 Mio particles selected using the resolution-enhancement 
method of Section [37T| while for halo models with f3{r) as in equa- 
tion l ll4t 2 Mio equal-mass particles were used. In order to ensure 
a careful modelling of the satellite, it was integrated with a shorter 
time step than most halo particles and the mutual forces with the 
halo particles were approximated with a much reduced opening 
angle. One simulation over 1000 time units (16000 block steps 
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Figure 2. Time evolution of the satellite orbits and halo Lagrange radii for four bound (left) and four parabolic satellite orbits (right) decaying in a halo with 
isotropic velocity distribution - note the different time scales. The top and middle panels show the evolution of the satellite distance to the halo centre and 
the satellite orbital energy Sj, respectively. Also shown in the middle panel is the halo's central potential depth il>[,(0) (lower curves). The thin vertical lines 
indicate the time at which the difference Ss - "thCO) is reduced by a factor 50. The bottom panel shows the evolution of the halo Lagrange radii (w.rt. the halo 
centre) containing 0.08%, 0.16%, etc. up to 20.48% of the halo mass. Within each set of orbits the initial orbital energy is the same, namely Es = for the 
parabolic orbits and Ss = -0.12 for the bound orbits, equivalent to that of the circular orbit at the halo half-mass radius (which is in fact the red orbit in the 
left panels). The purely radial orbits of both sets are plotted in cyan, while the colour sequence to red corresponds to ever less eccentric orbits, reaching e = 
(circular) for the family of bound orbits, corresponding to larger pericentric radii as indicated. 



or 256000 shortest time steps) took about 190 CPU hours (single 
CPU, N = I Mio), including some of the analysis. 

After each time step the position and velocity of the halo cen- 
tre was estimated from the position of the most bound particlefl 
Snapshots are stored at regular intervals and analysed in terms of 
their radial profiles using two different approaches. The first em- 
ploys simple averages over radial shells, assuming the centre to 
coincide with the one found from the most bound particles. The 
second estimates for each particle the density as a kernel estimate 
from its 32 nearest neighbours and then computes the centre, ra- 
dius, and other properties from density bins. This latter method is 
more robust in case the configuration is non-spherical, either be- 
cause of flattening or because of a spatial or velocity offset of the 
inner w.r.t. outer regions. 

For all except one halo model, simulations in isolation main- 
tained the original density profile over 2000 time units, except for 
the very inner parts, where artificial two-body relaxation and force 
softening result in a slight expansion. More quantitatively, the La- 

' Let denote a halo particle's specific potential energy due to all other 
halo particles (but not the satellite), then we first find the particle with 
the smallest (most negative) <f> and its K = 256 spatially nearest neighbours. 
Then we obtain the centre position and velocity as weighted average from 
the K/4 most bound particles (with largest |i^|) using the weights |0a:/4-i^,P. 



grange radii containing < 1 x 10"^ of the total mass expand noti- 
cably (see bottom panels of Fig. |2j, which is more pronounced in 
halo models with radially biased velocity distributions. The excep- 
tion is the halo model with initial /?(r) following equation dl4b with 
ySoo = 1 . This model turned out to be unstable to the radial-orbit in- 
stability and spontaneously re-arranges into a triaxial configuration 
within ~ 200 time units. In the course of this process the radial mass 
distribution is also slightly altered even without infalling satellite. 
In order to minimize the effects of these problems when interpret- 
ing our results, we compare each simulation with infalling satellite 
to a control simulation without satellite but identical initial halo. 
Table[T]shows a numerical summary of our results. 



4 ORBITAL DECAY 

In Fig. [2l we plot the time evolution of the satellite orbital radius 
and energy (top and middle panels) and the halo Lagrange radii 
(bottom panels) for both sets of satellite orbits in a dark-matter halo 
with isotropic velocities. In all these simulations the satellite mass 
and size are at their default values of nis = 0.01 and rs = 0.03. 
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4.1 Initial orbital decay 

We like to start our discussion by comparing tlie circular orbit start- 
ing at r = 4 (red in the lef t panels), which is similar to that used by 
IJardel & Sellwood ( l2009l) . and a plunging, almost radial, orbit with 
vanishing orbital energy (blue in the right panels of Fig.|2}. While 
both orbits decay to the centre, their evolution is clearly different. 
The circular orbit decays first slowly then faster, whereby remain- 
ing near-circular. The rate of decay, as indicated by the energy loss, 
is continuously increasing as the orbit comes closer to the centre. 
This is qualitatively consistent with Chandrasekhar's formula Q, 
since the halo density is steeply increasing inwards. 

The plunging parabolic orbit (blue in the right panels) on the 
other hand, suffers noticeable energy loss only near peri-centre, un- 
dergoing a stepped rather than steady decay of orbital energy. At the 
first peri-centric passage, the satellite loses enough orbital energy 
to find itself on a bound but highly eccentric orbit which returns to 
the inner parts of the halo within about 4 halo half-mass dynamical 
times. The energy loss at the second pericentric passage is larger, 
resulting in a quick orbital decay thereafter. The time scale (note 
the different the time axes in the left and right panels of Fig. |2j for 
the decay of this orbit is only about twice as long as that of the cir- 
cular orbit starting much closer (a circular orbit starting at r, does 
not decay within 1000 time units). Thus, even though more energy 
has to be lost for this orbit, the dynamical friction at peri-centre 
is so strong that the decay is still quite fast, even though the orbit 
spends about half its time well outside the halo virial radius. 

During the first few periods, the peri-centric radius for this 
plunging orbit is hardly decaying. This is expected if dynami- 
cal friction causes a near-instantaneous deceleration at peri-centre, 
which transfers the satellite to an orbit with the same pericentric 
radius. The apparent decay of the peri-centric radius after two pe- 
riods is not because dynamical friction away from pericentre be- 
comes significant, as the closest approach to the origin (as opposed 
to the halo centre) is in fact increasing, but because the centre of 
the halo has moved, as a consequence of the satellite interaction. 

For all orbits, we identify as the orbital decay or infall time 
finfaii the time at which the difference - 't>i,(0) first obtains 2% 
of its original value. There is a systematic trend of the infall time 
to decrease for higher eccentricities (smaller rp„i within each set of 
simulations), which is easily understood by the fact that dynamical 
friction is stronger for the higher density at smaller radii . 



4.2 Late orbital decay 

Interestingly, after the initial decay, when the orbital energy has 
almost reached its final value, the satellite and the halo centre are 
still orbiting each other at a distance of initially r ~ 0. 1 . In other 
words, the orbital decay was incomplete and has stalled. Unlike the 
situation investigated by previous authors, this stalling occurs at a 
radius which contains much less dark mass than the satellite mass 
and it seems more apt to say the halo centre orbits the satellite. 
This small orbit remains inert only for a short while (when 
\Xs - remains constant or even increases at ~ 0.1, see also 
Fig.ll3t. which is longer for high eccentricities of the initial orbit, 
while for the initially circular orbit (rpe,] = 4) there is hardly any 
such stalling. After this brief pause, the orbit decays very nearly 
exponentially with radial e-folding time of ~ 90 time units for the 
bound orbits (left panels of Fig. O, which is discemable for over 
two e-foldings until the orbital radius reaches the noise level. In 
case of initially parabolic satellite orbits (right), the stalled orbits 



are more eccentric and their decay more erratic (but qualitatively 
similar) with a longer decay time (e-folding time of ~ 160). 

Remarkably, the secondary decay times are very similar be- 
tween orbits within either set of simulation, but differ between 
them. This implies that the eccentricity of the initial and also the 
stalled orbit are not affecting the process responsible for this phe- 
nomenon. The differences in secondary decay times may be caused 
by differences, induced by the orbital decay, in the structure of the 
central halo, as discussed below. 

This secondaiy decay is not associated with any significant 
orbital energy loss and has no measurable effect on the final halo 
density profile (except perhaps for the innermost 0. 1 % of the halo 
mass), which is the main interest of our study. However, it is cer- 
tainly an interesting stellar dynamical phenomenon deserving fur- 
ther investigation. 



5 EFFECT ON THE HALO 

We now discuss the changes induced by the satellite's orbital decay 
in the dark-matter halo, concentrating mostly on the spatial distri- 
bution, while the changes to the velocity structure are mostly dis- 
cussed in the next section. 

5.1 Halo expansion 

The middle panels of Fig. [2] also show the evolution of the central 
halo potential <I'h(0) (lower curves; not including the contribution 
of the satellite), which at late time coincides with the satellite's or- 
bital energy. In all simulations, the final central potential depth of 
the halo is considerably shallower than initially, indicating a signif- 
icant reduction of the central dark-matter concentration. A compar- 
ison between the various simulations also shows that this reduction 
is more pronounced for parabolic than for bound satellite orbits, 
exactly as our analytic arguments of Section|2]predicted, as well as 
for more eccentric orbits. 

In the bottom panels of Fig. |2l we plot the time evolution of 
the halo Lagrange radii, while Fig. |3] shows for the same simula- 
tions the halo density p before and after the satellite infall, as well 
as the change ([9} in the cumulative halo mass profile and in the 
pseudo phase-space density p/cr^. From the time evolution of the 
Lagrange radii, we see that the expansion of the inner halo occurs 
rather suddenly at the time when the satellite settles in the cor^E 
in particular for the initially circular orbit (red in the left panels of 
Fig-IU. For more eccentric orbits, the early peri-centre passages re- 
sult in some minor expansion of the halo at radii comparable to the 
peri-centre radius, but hardly affect the innermost halo. 

The mass AM(r) removed from inside radius r (compared to 
the control simulatior^) has different radial profiles for the vari- 
ous simulations. Its amplitude is generally larger after the decay 
of a parabolic than a bound orbit, which is due to the larger or- 
bital energy of the former leading to stronger heating of the dark- 
matter particles. In fact, the maximum mass excavated AMn,ax is 
about twice the satellite mass if the latter is decaying on a plung- 
ing parabolic orbit, consistent with our models of Section |2l while 

^ The initial slow rise of the innermost Lagrange radii is entirely due to 
artificial two-body relaxation and not present in simulations with ten times 
the number of halo particles, see also the last paragraph of 33.31 
^ For the simulations, AM(r) is always measured this way, rather than 
against the initial model (as in equation |9), in order to account for halo 
evolution in absence of any satellite. 
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Figure 3. Radial profiles of the halo's pseudo phase-space density p/o"^ (top), density p (middle), and the change in cumulative halo mass AM (bottom) at 
f = 1000 and 2000 for simulations of the decay of bound (left) and parabolic (right) satellite orbits, respectively. The simulations and colour coding are the 
same as in Fig. |2] For each model, the vertical lines in the middle and bottom panels indicate the locations of rso% and r,nax. respectively, while the arrows 
in the top panels indicate the radius rM=nh within which the dark mass equals the satellite mass. The dotted curves correspond to the situation before satellite 
infall (f = 0), while the solid black curves represent the control simulation (no satellite). In the bottom panels, the black curves give M(r), not AM(r). The thin 
black lines in the middle panels are power-laws r^'^ ^(left) and r"" '' (right). 



for a circular orb it decaying from r = 4, c orresponding to the sit- 
uation studied bv'jardel & SellwoodI ( |2009|) . AM^ax < m^- This dif- 
ference between AM„,;,x obtained for orbits with Ss = (parabolic) 
Bs = -0.12 (bound) is even more pronounced than for the analytical 
models of Fig.[T] 

Within each set of orbits AM(r) is most peaked for the purely 
radial orbit, while it widens for lower eccentricities. This is pre- 
sumably because a satellite on a highly eccentric loses its energy 
to a narrow range of dark-matter particles close to its peri- centre, 
while less eccentric orbits lead to a distribution over a wider range. 



5.1.1 Central density reduction 

The density profiles in the middle panels of Fig. |3] confirm ear- 
lier results that the infall of a baryonic clump can considerably 
weaken the central dark-matter cusp. The final halo mass profiles 
(not shown) are noticeably perturbed interior to the radius r^ax (in- 
dicated by the vertical lines in the bottom panels of Fig.|3}, at which 
by definition the final halo density equals that of the control simu- 
lation without satellite. In order to quantify further the properties of 
the inner halo density profiles, we also calculate the radius rso% (in- 
dicated by the vertical lines in the middle panels of Fig. |3} interior 
to which the dark mass mass is reduced to 50% compared to the 
control simulation. In general, rso% is smaller than and marks 
the radius at which significant changes to the halo density profile 
are evident. 

There is a clear trend for more eccentric orbits to more 
strongly reduce the central density and hence result in a larger ra- 
dial range over which the density has been significantly reduced: 



r5o% is larger for smaller Tperf. The opposite is true for the radius 
''max- This contrasting behaviour of rso% and r,nax can be understood 
in terms of the different shapes of the excavated mass distribution 
AM(r), which becomes broader with increasing rperi, presumably 
because the energy input is spread over a larger radial range . 

Such an anti-correlation between r5o% and r,^ax for simulations 
with the same initial satellite orbital energy Es was also present 
in the analytic models in Fig. [T] For those analytic models the 
elfect was generated by assuming different central density slopes 
7o = -dlnp/dlnr|r=o- However, the central density slope of the fi- 
nal A'-body models (middle panel in Fig.[3]l are remarkably similar 
at around 70 = 0.4 - 0.5 at radii r < rso%- We should stress, how- 
ever, that the inner regions of these halo models are gravitationally 
dominated by the sunken satellite, and hence cannot be sensibly 
compared to the dark-matter haloes of galaxies. 



5.1.2 Central phase-space density reduction 

We also show in Fig.[3]the radii r/n^,,,^ at which the enclosed dark 
mass equals that of the satellite (arrows in the top panels). Evi- 
dently, these radii do not differ much from r^o%, implying that the 
innermost final dark-matter profiles are not self-gravitating: their 
gravity is dominated by the accreted satellite. Since the satellite 
size r^ = 0.03 is about ten times smaller than the radius inside of 
which it dominates the dynamics, it is effectively acting like a cen- 
tral point mass. A tracer population orbiting a point mass has a 
cusp p oc r"^'- if its phase-space density is constant. However, the 
simulations have much shallower central cusps, suggesting that the 
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Figure 4. Radial profiles of the intermediate-to-major {h/a, dashed) and minor-to-major (c/ci, solid) axis ratios, as well as the direction of the minor axis, 
provided as the .v, y, and z components of the unit vector along the minor axis (plotted where c/a < 0.95) for the final halo in the same simulations as in Figs.|2] 
and[3] The vertical lines in the top panel indicate for each model the location of rmax- 



actual dark-matter phase-space density has a central depression, de- 
creasing towards the highest binding energies. 

This is also borne out by the plots of the pseudo phase-space 
density p/<r^ in the top panels of Fig. [3] while initially and in the 
control simulations (dotted and solid black curves) these follow a 
pur e power law (this is how these models are actually constructed, 
see lDehnen & McLa ughlin 2005), the final profiles show a strong 
central depression, suggesting a considerable reduction of dark- 
matter phase-space density. The corresponding reduction in spa- 
tial density is much weaker, because the available phase space is 
increased above that from the self-gravitating cusp by the deeper 
central potential due to the accreted satellite. 

The central velocity dispersion initially decreases towards 
r = as cr^ oc ,-ni'n(r.2-r) fQj. ^ self-gravitating p oc r"'' cusp. In all 
simulations, the final halo has cr(r) increasing towards smaller radii 
(not shown, but evident from the strong central depression in p/cr^ 
compared to that in p), inside the radius rM=mf , as is required for any 
equilibrium system dominated by a central mass concentration. 

5.2 Halo shape 

In all our simulations, the dark matter halo is initially spherical. 
The top panels of Fig. |4] show the run of the final halo's principal 
axis ratios for the same simulations as in Figs. |2] and |3] At r < 
1.5 > Tniax the final dark matter distribution is near-oblate in all 
cases, becoming flatter towards the centre reaching c/a ~ 0.5 at the 
smallest measurable radius. The bottom panels of the same figure 
show the direction cosines of the minor axis. Note that the satellite 
initially orbits in the xy plane, starting at v = and x = n. For all 
but the purely radial orbits the halo minor axis is perpendicular to 
the initial orbital plane of the satellite. This is easily understood to 
originate from the transfer of orbital angular momentum from the 
satellite to inner halo particles during peri-centric passages, which 
presumably also generates the oblate inner halo shape. 

For the purely radial satellite orbits (cyan), which start off 
with zero orbital angular momentum, the final halo minor axis does 
not align with the satellite orbit, but appears to point in some ran- 
dom direction (though it is self-aligned). This behaviour is counter- 
intuitive as the initial models for these simulations are completely 
symmetric w.r.t. the infall axis. However, such a break of symme- 
try is the natural behaviour of radial orbits in the gravitational po- 



tentials of a density cusp. This is best explained by considering 
the limit of vanishing peri-centre radius for the change in azimuth 
A(p occurring over one radial period from apo-centre to apo-centre. 
For a harmonic potential, corresponding to 7 = 0, the radial orbit 
just passes straight through {A(p = n), while in the potential gener- 
ated by a point mass, corresponding to y = 3, the radial orbit is re- 
flected (A0 = 2n). For mass distributions with intermediate values 
of 7, such as dark-matter haloes, Aip is between these two extremes, 
and the synmietry of the initial orbit is brokerQ 

This deflection of the satellite orbit is compensated by an 
equal and opposite momentum to the inner halo, such that the sub- 
sequent relative orbit of the two is no longer radial. The further 
evolution follows the same pattern as for the non-radial orbits: the 
halo flattens perpendicular to the angular momentum axis of this 
orbit. Note that for both radial orbit simulations shown in Fig. |4] 
the minor axis is near-perpendicular to the original infall direction 
(x-axis). This is explained by the fact that A(p for a shallow cusp 
is only slightly larger than the value n for the harmonic potential, 
implying that the orbit is only weakly deflected from its original 
infall trajectory, which is indeed what we see in our simulations. 

There is also some hint of triaxiality at small radii, in partic- 
ular for the radial bound orbit. This is somewhat surprising as it 
mostly occurs within the region where the satellite dominates the 
enclosed mass and hence the potential is near-spherical. We suspect 
that this very central triaxiality is generated during the secondary 
orbital decay if the small decaying orbit is eccentric. 



6 EFFECT OF HALO VELOCITY ANISOTROPY 

So far, we have discussed the results from eight simulations, which 
differed only in the satellite orbit, but had the same initial halo 

* Strictly speaking, the radial orbit with L = is never deflected (Ai^ = n) as 
no transverse forces act on it. However, in the limit L — > one gets Aip > n 
(except for the harmonic potential), i.e. is discontinuous at L = or, 
equivalently, r^i.^\ = 0. In our simulations, the halo potential is modelled 
from the softened particle potentials and deviates from the power-law form 
at r < e = 0.005, where it becomes harmonic and the discontinuity at rpe^ = 
is removed such that A(p x n for rpe,i < e. However, even the simulations 
with initially purely radial satellite orbits have actual rpe,i > e at first pas- 
sage (see Fig.|2), such that the simulated satellite orbit is actually deflected. 



10 David R. Cole, Walter Dehnen, Mark I. Wilkinson 



0.01 



. 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
; /3,=-0.43 




1 : 


: /3i= 




- 


- /3, = 0.3 






P{r) ^ 




\ : 






\ 

; \ 

\ \ - 
\ - 








- i 

:h 






■ ll 
■.ill 

It 

f , , , , 1 , , , , 1 , , , , 1 




1 



-0.5 



-0.4 



-0.3 -0.2 
E 



-0.1 



Figure 5. Differential energy distributions for the four different halo models 
considered, which only differ in their velocity anisotropy profiles. 
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Figure 6. Radial profiles of the final (at t = 2000) halo velocity anisotropy 
parameter/? for 4x3 simulations: four parabolic satellite orbits (as shown 
in the right panels of Figs. |2]to|4] with the same colour coding) decaying 
in three halo models with different constant initial anisotropy (line style as 
indicated). The black and green curves are the radial anisotropy profiles 
(obtained in the same way from W-body data) of the initial conditions and 
the control simulations, respectively. The satellite mass and size are m, = 
0.01 and rj = 0.03 (as for all simulations presented sofar). 



model with velocity isotropy. We are now investigating initially 
anisotropic dark-matter velocity distributions. Since parabolic or- 
bits are more realistic for infalling clumps, we restrict our discus- 
sion to the four satellite orbits with initially vanishing orbital en- 
ergy Bs = and various eccentricities. Most of our general conclu- 
sion are, however, at least qualitatively also valid for bound orbits. 

In addition to the halo model with velocity isotropy, used 
in Sections |4] and [5] we consider three initially anisotropic mod- 
els. Two of these have constant anisotropy parameter of y6 = 0.3 
and j3= - 0.43, respectively, which corresponds to the same 
level of anisotropy in the sense of |lno"r/crt| (with the tangen- 
tial velocity dispersion cr^ = [cr^ + cr^]/2). The third anisotropic 
halo model is motivated by simulations of galaxy halo formation, 
which generally predict that the velocity distribution of the dark 
matter wi thin haloes is outwardly increasing radially anisotropic 
( [Hansen & Mo ore 2006). We use the anisotropy profile of equation 
il4\ withjSoo = 1. This halo model is quite different from the others, 
as it undergoes a radial orbit instability, and we discuss it in some 
more detail in section|7] 



6.1 Halo phase-space structure and vulnerability 

Before we discuss the simulations, let us consider an important dif- 
ference between these four models. Haloes with more radially bi- 
ased velocities have more mass on eccentric orbits, which spend 
most of their time near apo-centre, but contribute significantly to 
the density near their respective peri-centres, in particular if the 
density increases more slowly than r^^ (because of geometrical ef- 
fects) as it does in the inner regions of dark-matter haloes. Conse- 
quently, much of the dark matter in the innermost region of a halo 
with radial velocity anisotropy is on orbits which spend most of 
their time outside the innermost halo and have lower binding en- 
ergies than the local circular orbits. Thus the more radially biased 
the velocities, the less mass is at high binding energies, as demon- 
strated in Fig. |5] which shows the differential energy distributions 
dM/dE for our four halo models (computed from W-body data; for 
a cleaner plot of the same effe ct but for different models see Fig. 4.5 
of lBinnev & Tremainell20o3) . 

The relative lack of highly bound orbits in haloes with ra- 



dial velocity anisotropy has immediate consequences for the re- 
sponsiveness and hence vulnerability of the central regions to 
perturbations, such as a n infalling massive satellite (see also 
iBinnev & Tremainell20oi p. 299). Highly bound orbits are con- 
fined to the central regions, which makes them relatively inert to 
external perturbations, in particular if their orbital period is short 
compared to the time scale of the perturbation (adiabatic invari- 
ance). The eccentric orbits in haloes with radial velocity anisotropy, 
on the other hand, have longer periods and hence are not adiabat- 
ically protected in the same way. Moreover, the infalling satellite 
can relatively easily perturb dark-matter particles near the apo- 
centres of such orbits, increasing their angular momenta and hence 
peri-centric radius, and thereby reducing the central density of the 
halo. 

These arguments also suggest that the difference in vulnera- 
bility between haloes with isotropic and radially anisotropic veloc- 
ities is more pronounced for perturbation by a satellite passing at 
larger peri-centric radius, which affects eccentric orbits contribut- 
ing to the centre, but hardly the innermost orbits. A satellite falling 
in on a purely radial orbit, on the other hand, will affect dark-matter 
orbits at all radii, regardless of halo anisotropy. 

6.2 Change in velocity anisotropy 

In Fig.[6l we plot the radial y6 profiles of the final halo in simulations 
of the decay of the four different parabolic satellite orbits (colour 
coded as in Figs.[2]to|4]l in each of the three different halo models 
with constant initial p. 

In case of initial velocity isotropy the final velocity distribu- 
tion is isotropic too, while in all of the initially anisotropic cases, 
the halo velocity distribution evolves towards isotropy in the inner 
regions. This evolution is partly driven by (artificial) two-body re- 
laxation, as evident from the control simulations (green in Fig.|6]l, 
which results in velocity isotropy at r < 0. 1. However, for the sim- 
ulations with decaying satellite orbit, this isotropisation occurs to 
larger radii, which is essentially independent of the orbital eccen- 
tricity of the decaying orbit, and is only a function of the initial halo 
anisotropy. 

This radius is larger by about a factor 3 for the models with 
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Figure 7. Dependence of orbital decay time fi„faii and several properties of 
the final halo (the bottom panel plots the radius inside of which the minor- 
to-major axis ratio c/a < 0.8) on initial halo velocity anisotropy (hne style 
as in Fig. [6) and initial peri-centre for parabolic orbits with satellite of mass 
fjis = 0.01 and size = 0.03. Symbol colours match those in Figs.|2]to|4] 



initially constant tangential velocity anisotropy than for the mod- 
els with initially constant radial velocity anisotropy. This is some- 
what surprising in view of the discussion in the previous subsection. 
One speculation is that a system with tangential velocity anisotropy 
is not well-mixed in the sense of iTremaine. Henon & Lvnden-Belll 
( Il986h and Dehnen ( 2005), such that violent relaxation, induced by 
the perturbation, promotes evolution towards isotropy, while per- 
haps the opposite is true for radial velocity anisotropy. 

It is also interesting to note that the evolution of the anisotropy 
profile is complete by the time of the first qualitative change in 
the orbit (i.e. essentially at / = /infan). Thereafter, the final secondary 
decay of the satellite orbit takes place without further modification 
of the halo velocity distribution. 



6.3 Effect on orbital decay and final tialo 

Instead of showing detailed figures, similar to Figs.|2]to|4] of the 
time evolution and radial profiles of density, axis ratios, etc. for the 
final halo of all the 4 x 4 simulations, we summarise the effects 
of different initial halo anisotropy in Fig.|7] which shows, for each 
simulation, the infall time ?i„faii as well as several properties of the 
final halo (at / = 2000). Apart from AMy„.,,^/mg, r^^^^, and r^o%, al- 
ready used in previous sections, we also plot the maximum value 




1500 



Figure 8. Like the right panels of Fig.|2] but for an initial halo with cos- 
mologically motivated velocity anisotropy profile. 



for the pseudo phase-space density p/cr^ and the radius inside of 
which the minor-to-major axis ratio c/a < 0.8. 

There is a systematic trend of shorter infall times for initially 
more radial dark-matter velocity anisotropy, exactly as expected 
from the arguments of Section 16.11 the inner parts of haloes with 
radially biased velocities are more responsive to and vulnerable by 
the infalling satellite. This is also the reason for the behaviour of 
^50% and max{p/cr^), which demonstrate that haloes with initially 
more radial velocity anisotropy suffer the strongest reduction in 
their central density and phase-space density. As argued in Sec- 
tion l6.1l the effects of initial halo velocity anisotropy are most pro- 
nounced for satellite orbits with large rpj-^ and least for a purely 
radial satellite orbit. 

Again, the radius r„,a,K at which AM(r) peaks anti-correlates 
with the radius rso%, indicating that the distribution AM(r] becomes 
more peaked (not shown) for more radially biased velocities and 
that this peak occurs at smaller radii. Both of these are natural con- 
sequences of the difference in orbital structure as outlined above. 

There is also a systematic effect on the halo shape. For each 
model with initially constant /?, the bottom panel of Fig.|7]plots the 
radius inside of which the minor- to-major axis ratio c/a < 0.8. The 
halo with tangential velocity anisotropy becomes flattened in a re- 
gion comparable to that of the halo with isotropic velocities. For the 
halo with yS = 0.3, the flattening is more pronounced, reaching out 
to about twice as far, and the shape is near-oblate (not shown) with 
the same characteristic as for the halo with initial velocity isotropy, 
discussed in Section \5l2\ This is again expected, as the satellite's 
orbital angular momentum absorbed by the halo is relatively more 
important for dark-matter on low-angular- momentum orbits, as for 
radial velocity anisotropy. 
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Figure 9. Final halo radial profiles after the infall of a satellite on an ini- 
tially parabolic orbit into a halo with initially radially increasing velocity 
anisotropy /?(r) as in equation [14] with /?„ = 1- The plots are equivalent to 
Figs.[3]and|4](the thin line in the density plot is r"""^), except that in the bot- 
tom three panels we plot not only the minor-axis direction cosines (solid), 
but also the major axis direction cosines (dashed). For this halo model the 
control simulation (black) also undergoes some evolution driven by a radial- 
orbit instabihty, most evident from the changes in anisotropy and shape, see 
text for details. 



7 HALO WITH COSMOLOGICALLY MOTIVATED 
VELOCITY ANISOTROPY 

The situation with a outwardly increasing radial velocity anisotropy 
as in equation jl4b is typical for dark-matter haloes emerging from 
simul ations of large-scale structure formation (Hansen & Moore 
l2006h . As already mentioned in section 13.31 our spherical halo 
model with such an anisotropy profile undergoes a radial-orbit in- 
stability and quickly settles into a prolate configuration. While this 
instability somewhat complicates the inteipretation of any results, 



the corresponding simulations are presumably most realistic what 
concerns the halo structure. 



7.1 Orbital Decay 

For the four simulations of initially parabolic satellite orbits in this 
halo model, we plot in Fig.[8]the orbital decay and time evolution 
of the halo Lagrange radii, similar to the right panels of Fig.|2] First 
note that the control simulation (black) undergoes an initial ex- 
pansion within the first ~ 200 time units. This expansion is largely 
driven by the violent relaxation during the re-arrangement to a pro- 
late configuration (due to the radial-orbit instability), while the ex- 
pansion seen in the other control simulations (see Fig.|2) was solely 
driven by (artificial) two-body relaxation. 

The orbital decay is faster than for any other halo model con- 
sidered (see top panel of Fig. |7}, because of the more efficient 
transfer of satellite orbital energy and angular momentum to the 
halo particles (as outlined in Section lOl during the very first peri- 
centric passages, when the halo density is not yet diminished by 
the instability-driven expansion. The associated heating of the in- 
nermost halo results in a slightly faster expansion in the simulations 
with satellite than in the control simulation even at r < rperi . 

The same is true to a lesser degree for the halo with /? = 0.3 
(not shown), while for the halo with isotropic velocities (Fig. O, 
the Lagrange radii at r < rp„i were hardly affected by the early 
peri-centric passages. This different response of the innermost halo 
to the satellite's first peri-centric passage can be attributed to the 
predominance of eccentric dark-matter orbits in the innermost halo, 
which are perturbed by the passing satellite at their respective apo- 
centres, as outlined in Section|6.1| 



7.2 Effect on lialo structure 

Fig- m shows the radial profiles of p/cr^, p, AM, fj, c/a, b/a and 
the direction cosines of the minor and major axes (similar to the 
right panels of Figs. |3] and |4j for the four simulations with decay- 
ing satellite orbits (coloured), the initial model (dotted), and con- 
trol simulation (black). All these radial profiles after satellite orbit 
decay are quite similar, the strongest difference is 15% between 
the amplitudes of AM(r). This similarity is also expected from our 
discussion in Section [6T] the halo is so responsive that the exact 
satellite orbit does not matter too much. 

The pseudo phase-space density, p/cr^ (top panel), of the con- 
trol simulation is even larger than initially, although the density is 
reduced at r < 0.15. The reason for this counter-intuitive behaviour 
is that the isotropisation (evident from the runs of /3) has reduced the 
velocity dispersion cr in the inner parts. However, after the orbital 
decay of the satellite, p/tr' is substantially reduced and necessar- 
ily also the true phase-space density. The maximum pseudo phase- 
space density is a factor 2-3 smaller than for simulations with initial 
velocity isotropy, and the density a factor ~ 2. This represents the 
strongest central halo reduction in all models (for the default val- 
ues of satellite mass and size), which makes sense in view of the 
vulnerability due to the radial velocity anisotropy. 

The mass AM excavated compared to the control simulation 
is slightly less for the purely radial orbit (but slightly more for the 
orbit with j-peii = 1.3) than in case of an isotropic halo (Fig.|3}. 
However, such a comparison is not quite adequate in view of the 
different behaviour of the respective halo models in isolation, and 
the picture that stronger radial velocity anisotropy results in larger 
an effect on the halo remains valid. 
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Figure 10. The variation with satellite mass of finfaii. AM,„axi '■max, >'50%, 
and the maximum (over all radii) of p/cr' after the decay of a circular satel- 
lite orbit starting at r; = 4 in a halo with isotropic velocities (red) or after 
the decay of a parabolic satellite orbit with Tp^^ = 0.4 starting at r; = 10 in a 
halo with outwardly increasing radial velocity anisotropy (blue). The lines 
are power-laws with exponent as indicated. 



In all simulations is the final halo more isotropic than initially 
at r < 3, but still retains significant radial anisotropy. In fact, the j3 
profiles are identical to that of the control simulation, presumably 
because this profile corresponds to a well-mixed state, attained after 
the initial violent relaxation phase, and hardly alTected by any fur- 
ther relaxation due to satellite interaction. Such radially anisotropic 
velocity distributions are likely the inevitable property of triaxial 
and prolate equilibria, because of the p redominance o f low-angular- 
momentum orbits, such as box orbits ( lDehnenll2009l) . 



7.3 Effect on halo shape 

Obviously, the final halo shapes are completely different from those 
obtained in simulations with any other initial halo model. As a re- 
sult of the radial-orbit instability, the control simulation obtains 
a strongly flattened purely prolate shape with axis ratio < 0.6 at 
r < 2 (and some random orientation). This is in remarkable agree- 
ment with the fact that dark-matter haloes in gravity-only simula- 
tions o f large-scale structu re formation tend to have near-prolate 
shapes dWarren et al .11 19921) . After the satellite infall, however, the 
halo shape becomes less flattened and triaxial in the inner parts. 
The changes in halo shape compared to the control simulation ex- 
tend to about r ~ 5, much farther than the changes in halo density 
or velocity anisotropy. 

As the direction cosines indicate, the orientation of the princi- 
ple axes of the triaxial shape is constant with radius, as one expects 
for an equilibrium system. The minor axis is always perpendicular 
to the original satellite orbit, while the major axis is somewhere in 
the initial orbital plane (except for the purely radial orbit for which 
such a plane cannot be defined). 



Figure 11. Like Fig. 1101 except that satellite size is varied and nij = 0.01. 
kept constant. 



8 VARYING SATELLITE MASS AND SIZE 

All simulations presented sofar have the same satellite mass = 
0.01 and size = 0.03. We are now investigating effects of larger 
and smaller values for these parameters. 



8.1 Varying satellite mass 

Fig.[TO]shows the effect of varying by a factor of 3 up and down 
for a satellite initially on a circular orbit in a halo with isotropic 
velocities (red) and for a satellite on a parabolic with r^^i = 0.4 
in a halo with cosmologically motivated /3(r). As expected, more 
massive satellites reach the centre more rapidly, and cause more 
damage to the inner halo. In agreement with Fig.[T] r,nax, ''50%, and 
AM,^ax are largest for the most massive satellite, and decrease sys- 
tematically as the satellite mass is decreased. As the figure shows, 
the scaling of r„ax and rso% with satellite mass is close to scaling 
r oc m°-^ (solid lines) also found for our analytic models in Fig.[T] 

The most interesting result from this study of varying satellite 
mass is that low mass satellites are more efficient at displacing mass 
than high mass satellites: in agreement with the analytic models 
the relative excavated mass AM,^ax/»is increases towards smaller 
satellite masses. 

According to Chandrasekhar's dynamical friction formula ([T]l, 
the drag on the satellite is proportional to its mass and would 
therefore naively expect the infall time to scale inversely with satel- 
lite mass: t]„f^]] oc m^' . This is in fact the scaling indicated by the 
solid lines in the top panel of Fig. [TO] However, the actual infall 
time measured for our simulations differ from this expectation. For 
the simple situation of an initially circular satellite orbit in a stable 
halo with isotropic velocities (red triangles), the infall time scales 
more like t;„f.^n oc in^'^ '^. This difference is most likely caused by the 
re-adjustment of the halo during the satellite infall, an effect we ig- 
nored when deducing finfaii oc ot^t', which therefore only applies in 
the limit — > 0. The fact that this limit is not applicable even for 
OTs/Mhaio = 0.01 may seem surprising, but is not in view of the fact 
that only a small fraction of the total halo mass accounts for most 
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of its binding energy, whiich in turn enables such a feeble satellite 
to inflict significant damage to the halo. 

The orbital decay into the halo model with radially increas- 
ing f}{r) (blue squares) is complicated by the fact that the halo si- 
multaneously undergoes a radial-orbit instability. The slight cen- 
tral density reduction driven by this instability increases the in- 
fall time compared to the original halo within the first ~ 200 time 
units. In our simulations, this sets up a race between satellite in- 
fall and instability-driven halo evolution. A more massive satellite 
sinks more quickly than this evolution and t^^f.^w of the original halo 
applies, while for a low-mass satellite the longer infall time after 
the instability-driven halo evolution applies. Clearly, this race sce- 
nario is an artifact or our simulation setup, and we expect, based 
on the arguments given above, that in general /infaii scales slightly 
shallower than m^'. 

8.2 Varying satellite size 

In Fig.[TT]we vary the size of the satellite by a factor 3 up and 
down, while holding its mass fixed for a circular satellite orbit in 
an isotropic halo (red triangles) and a parabolic satellite orbit with 
''peri = 0.4 in a halo with cosmologically motivated p{r). There is 
a systematic trend that more compact satellites lose energy more 
quickly and hence exhibit more rapid orbital decay. This is rea- 
sonable as more compact satellites are more efficient at scattering 
background particles and thus lose energy to the halo more rapidly. 

The panels of Fig. [TT] show that the cumulative effect of the 
satellite on the halo is almost independent of its size - the values of 
AMi^ax/^si ''max 5 and maxjp/cr'j are virtually unchanged as the size 
of the satellite is changed by an order of magnitude, in particular 
for the parabolic orbit in the p{r) halo. This is not very suiprising 
in view of the analytical models of Section|2] the orbital energy of 
the satellite is essentially independent of (as long as ^ r2). 

However, for the circular satellite orbit decaying in a yS = 
halo, there is a systematic trend of larger r5o% with larger and 
the halo profile is more significantly flattened by a more extended 
satellite. This is also seen in the evolution of the Lagrange radii and 
is reasonable since the more extended the satellite, the less energy 
it loses at large radii and therefore has more energy to impart to 
the innermost halo. More extended satellites on near-circular orbits 
carry more energy to the inner regions and hence perturb the halo 
profile to a greater degree (for highly eccentric orbits this picture 
does not apply, as they affect the innermost halo already at their 
first peri-centric passage). 

The evolution produced by satellites on bound orbits in 
isotropic haloes is qualitatively the same as that for parabolic or- 
bits, although the magnitude of the effects is reduced. In particular, 
the infall time differs by only 30% between the largest and smallest 
satellites, compared to the factor of three for the parabolic satellite 
orbit in the halo with increasing pir). This again may be caused by 
the increased vulnerability and responsiveness of haloes with radial 
velocity anisotropy, as discussed in Section|6T| 



9 THE EFFECT OF SATELLITE REMOVAL 

In our modelling so far we have ignored the additional halo expan- 
sion following the possible loss of baryons in a feedback-driven 
galactic wind. At the final time in the models, the gravitational po- 
tential inside ~ 0.1 r2 is dominated by the satellite. Therefore the 
central potential is still quite deep, even though for most of our 
models the satellite mass is somewhat smaller than the removed 




r 



Figure 12. As Fig. |3] but also including the removal of the satellite. The 
blue curves refer to the same model as the blue curves in Fig. |3] The ma- 
genta curves are obtained from this model after instant removal of the satel- 
lite; while for the red model the satelite has been slowly removed. 

dark mass. This implies that the bound phase-space volume avail- 
able for the dark matter at the centre is almost unchanged. 

In order to study the effect of baryon outflow on the dark- 
matter, we extend one of our models by removing the satellite. The 
model concerned is the one with an initially parabolic satellite or- 
bit with Tpeii = 0.4 falling into a halo with isotropic velocities {blue 
curves in the right panels of Figs. |2]and|3). In Fig.[T2] the result- 
ing radial profiles for pier", p, and AM are shown for the situation 
(i) after satellite accretion (same as the curves in Fig. |3] blue), (ii) 
further 1000 time units after instant satellite removal (correspond- 
ing to a fast wind; magenta), and (iii) after a slow satellite removal 
(over 250 times units followed by another 250 time units without 
satellite; red). Not surprisingly, the satellite removal has a dramatic 
effect on the central halo density: the total mass excavated has al- 
most doubled compared to the situation prior to satellite removal 
and the density has become clearly cored with a central density 
~ 10 times smaller than the initial density at the radius where ini- 
tially = M{r). 

The most intriguing plot, however, is that of the pseudo phase- 
space density p/cr^ {top panel of Fig. [72l >. While the re-distribution 
of material to larger radii has shifted the iimer profile somewhat, 
the overall distribution and the maximum value for pla^ remained 
unaffected by the satellite removal. For the slow-wind model this 
is, of course, expected, since the (true coarse-grained) phase-space 
density is conserved. Therefore, the additional reduction of p and 
AM in this case is entirely owed to the reduction in available phase- 
space effected by the removal of the satellite potential. Of course, 
this is a non-linear process because the reduction in central dark- 
matter density itself leads to a further decrease in the potential 
depth. This explains why the additional mass removed by the satel- 
lite departure (the difference between the red/magenta and blue 
curves) is quite substantial, exceeding the satellite mass. 
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Figure 13. Secondary orbital decay: time evolution of the distance between satellite and halo centre for the 300 time units after the infall time (vertical line) 
for the bound satelhte orbits in halo models with various velocity anisotropy as indicated (the curves for different halo models are offset by 2dex from each 
other). The horizontal dashed lines indicate the noise level measured at much later times. 

\Xs - Xsl between halo centre and satellite after the initial infall of 
the satellite. In all these cases the secondary orbital decay occurs. 
The secondary orbit may be eccentric (as for r^^^i = in Fig.ll3t or 
near-circular, when the decay appears more regular. In the case of 
parabolic satellite orbits (not shown), the secondary decay does not 
always occur or is much more noisy and less regular. 

We are not sure about the cause and dynamics of this phe- 
nomenon, and whether it is a numerical artifact or not. It may be re- 
lated to the form of the softened satellite potential, though the sec- 
ondary orbital amplitude decays from ~ 0.2 > to ~ 0.007 < r,, 
such that the harmonic inner parts of the softened potential become 
relevant only in the late stages. 



Somewhat unexpectedly, the fast-wind model is only slightly 
more efficient at removing dark-matter. This strongly suggests that 
even a fast galactic wind is nearly adiabatic, i.e. has little effect on 
the (true coarse-grained) phase-space density (which is unaffected 
by a slow wind). 



10 SUMMARY 

We have modelled, using A'-body simulations, the decay of satellite 
orbits in spherical dark-matter halo models with various degrees of 
velocity anisotropy. 



10.1 Orbital evolution 

The evolution of the satellite orbits has several phases. First there is 
a steady decline in the amplitude of the orbit with almost constant 
peri-centric radius, driven by dynamical friction near peri-centre. 
This is followed by a period of rapid shrinkage of the peri-centre 
coinciding with significant expansion of the innermost halo, as ev- 
ident from the evolution of the Lagrange radii in Figs.|2]and|8] 

At this stage the evolution of the halo, driven by the transfer 
of energy and angular momentum from the decaying satellite orbit, 
is more or less complete. The time required for this process scales 
as finfaii 'WsT"' with satellite mass, somewhat shallower than the 
inverse scaling expected from Chandrasekhar's dynamical-friction 
formula. We also find that finfaii is significantly shorter for decay in 
haloes with radial velocity anisotropy, because the orbital structure 
of such haloes makes them more susceptible to perturbations. 

The orbital evolution shows some interesting and unexpected 
behaviour after fjnfaii. The most bound region of the halo and the 
satellite form a sort of binary, dominated in mass by the satellite, 
which is modelled as a single softened particle with size (soften- 
ing length) =0.03 (our unit of length is the halo scale radius). 
Initially, this secondary orbit has amplitude x 0.1 > rs, but after 
a brief period of apparent growth decays very nearly exponentially 
with e-folding time of 90 time units. 

This is clearly visible in Fig. [13] which plots for all 12 simula- 
tions of bound satellite orbits the late time evolution of the distance 



10.2 Effect on Halo 

Even though the assumed satellite mass was only 1% of the total 
halo mass, the damage done to the halo is significant: the inner 
parts of the halo are substantially reduced in density and phase- 
space density (as indicated by the behaviour of pjcr'), when the 
satellite has displaced up to twice its own mass from the innermost 
halo. The satellite also affects the velocity structure of the inner 
halo making it more isotropic. Finally, the halo shape becomes tri- 
axial in its inner parts. 

We find that the efficiency with which the satellite can affect 
the halo, e.g. the amount of central density reduction, is increased 
for haloes with radial velocity anisotropy. This is understandable in 
terms of the orbital structure of such haloes, making them more vul- 
nerable (see Section l6Tt . and relevant, as real dark-matter haloes 
are expected, from simulations of large-scale structure formation, 
to be radial anisotropic. 

One of our most interesting findings is that satellites of smaller 
mass are relatively more efficient in damaging the halo: the ratio 
AMmax/"'s of the maximum displaced halo mass over the satellite 
mass increases towards smaller satellite masses (though the or- 
bital decay time decreases, of course), as in fact predicted by the 
analytic energy arguments presented in Section |2] 
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11 DISCUSSION 

In this study we have ignored the details of baryonic physics and 
instead used the simple model of a compact baryonic clump falling 
into a dark-matter halo. Of all possible scenarios for the effect of 
baryonic physics on the structure of dark-matter haloes this is likely 
the most efficient. This process proved to be highly effective in al- 
tering the structure of the dark-matter halo in its inner parts, where 
the galaxy resides. A single clump of only 1% of the total halo mass 
can reduce the density at < 0.1 halo scale radii by more than an or- 
der of magnitude (including the effect of a galactic wind as demon- 
strated in section |9}, more than sufficient to explain the discrep- 
ancies between the rotation curves predicted (using gravity-only 
simulations) and observed for dark-ma tter dominated galaxies (e.g. 
ISimon etal ll2005l : Ide Blok et allzOOSl) . 

The total amount of baryons is one sixth of all matter 
iKomatsu et alj |201C|) . i.e. 0.2 times the dark mass, 20 times more 
than our model clump. Thus, if the baryonic heating of a dark- 
matter halo is only ~ 2% efficient on average, the damage done to 
the halo is equivalent to that in our models, because any additional 
addiabatic inffow and outflow of baryons has no lasting effect on 
the dark matter. 

This assumes, of course, that the baryonic sub-structures and 
clumps are able to heat the dark-matter particles which contribute 
to the central cusp. Our simple model of a compact clump, which 
sinks to the centre of the dark-matter halo, is certainly somewhat 
unrealistic as baryonic structures are susceptible to disruption by 
tidal forces before they reach the innermost parts of the halo. How- 
ever, owing to their dissipative nature baryons, unlike dark-matter, 
can form new sub-structures and clumps, which then continue to 
heat the dark-matter via dynamical friction. Moreover, as our sim- 
ulations have shown, radial velocity anisotropy, which is typical for 
dark-matter haloes, make their innermost parts more vulnerable to 
perturbations from the outside. This is because much of the matter 
in the cuspy region is on eccentric orbits spending most of their 
time at much larger radii where they are prone to dynamic heating 
(see also Section lOl . 

11.1 Dark-matter contraction vs. expansion 

A key process in the re-shaping of dark-matter haloes by non- 
gravitational baryonic physics is the transfer of energy via dynam- 
ical friction from baryonic sub-structures to the dark-matter parti- 
cle|3. The effect of this heating can be understood qualitatively by 
considering the Jeans equation of hydrostatic equilibrium 

Vipo-^) = -pVO. (17) 

The heating of the dark matter by the non-adiabatic baryon infall 
raises the central cr^, which, at fixed p, increases the pressure gradi- 
ent (left-hand side). At the same time, the arrival of the baryons in- 
creases the gravitational pull (right-hand side). If these two effects 
balance exactly, the dark-matter density p remains unaffected. If the 
heating dominates, as was the case in our maximally non-adiabatic 
simulations, then p must flatten to retain equilibrium. Conversely, 
if the gravitational pull dominates (especially for negligible heat- 
ing, i.e. 'adiabatic contraction'), then p has to steepen (and tr^ will 
increase adiabatically due to the compression). The exact balance 
depends on the details and most likely varies systematically with 

' While this process itself is, of course, purely gravitational, it is neglected 
in gravity-only simulations, which ignore the formation of baryonic sub- 
structures and unequivocally predict cuspy dark-matter haloes. 



galaxy type, size, and environment, explaining the possibility of 
conflicting results from si mulations which attempt to model baryon 
physics more directly (e.g. lRomano-Diaz et al.l2008l : |Pedrosa et al.1 
2009). 

An alternative way to look at this problem is to consider 
the effect on the dark-matter (coarse-grained) phase-space density, 
which arguably is dynamically more relevant than the density, be- 
cause it is conserved for adiabatic evolution. However, the process 
of baryon infall is inevitably non-adiabatic and reduces the dark- 
matter phase-space densitjtj. As long as the accreted baryons re- 
main at the centre, this may not necessarily result in a reduction of 
the dark-matter spatial density, because the additional gravitational 
potential of the newly arrived baryons increases the available bound 
velocity-space (phase-space at fixed position). 

Of course, a galactic wind changes the situation: the loss of 
some or all of the accreted baryons tips the balance towards a re- 
duction of the dark-matter density, as convincingly demonstrated 
by our models of satellite removal. In the Jeans-equation picture 
the wind removes the additional gravitational pull from the accreted 
baryons. In the phase-space interpretation, the wind reduces the 
bound velocity space, thus pressing the dark-matter phase-space 
fluid out of the centre, like toothpaste out of its tube. 

11.2 Application to dSph and GCs 

For a dwarf spheroidal galaxy our simplistic model of baryon in- 
fall may apply even more directly. Given that their present-day 
baryonic mass is comparable to the dark mass that needs to be 
rearranged in order to convert their haloes from cusped to cored, 
it seems plausible that for a reasonable star formation efficiency 
one could build the stellar component of the dSphs from a num- 
ber of star clusters that fall to the centre by dynamical friction 
and in so doing generate a macroscopic core in the halo dark mat- 
ter distribution of the dSph. For this to work, the clusters need 
to remain largely unscathed b y tides before they reach ~ 0.1r2. 
IPeiiarrubia. Walker & Gilmord (l2009b have investigated the tidal 
disruption of star clusters orbiting in a dSph, in particular the Sagit- 
tarius and Fornax satellites to the Milky Way, and used A'-body 
simulations to confirm the following criterion for the tidal radius r, 
of a globular cluster (GC) at radius r within a halo 

PGc(n) « 3ph(r). (19) 

In particular they looked at the resilience of the five GCs in Fornax 
and found that GCs which retain bound masses of greater than ap- 
proximately 95% of their total mass do not undergo tidal disruption 
by the halo. In our simulations the clump falls in to approximately 
r = 0.1 f2 at which point it has excavated the mass in the centre of 
the halo and reduce d the central densi ty. At a corresponding radius 
in the simulations of lPeriarrubia et all four of the GCs are still very 
stable against tidal disruption. The fifth one is only disrupted when 
it spends a large fraction of a Hubble time at a radius in the range 
corresponding to 0.1 - 0.2r2 of our model. 

In a similar way we can look a t the stability of GCs i n Low 
Surfa ce Brightness (LSB) galaxies. iKuzio de Narav et alj ( 1200 Sl 
I2OO6I) investigated the density profiles of a number of LSB galaxies 
and attempted to fit them using NFW and pseudo-isothermal halo 

* More precisely, it reduces the excess-mass function 

with f(x. v) the coarse-grained phase-space density toehnenlboOSb . 
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models. Their best-fit NFW models should give an upper limit to 
the density in the inner halo and therefore the one most likely able 
to disrupt an infalling object. Using these we find a range of val- 
ues for Pi,(r) at 0. 1 scale radii ranging from 0.012 to 0.080 MqPC"^^. 
These are lower than the correspondi ng value for Fornax whi ch is 
~ O.lSMopc"^ Based on the work of IPeiiarrubia et al.l (l2009h this 
implies that GCs would also be stable at 0.1^2 in LSB galaxies. 
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